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1. Introduction 



Abstract. We present reliable a-postcriori error estimates for /ip- adaptive finite element approxima- 
tions of eigenvalue/eigenvector problems. Starting from our earlier work on h adaptive finite element 
approximations we show a way to obtain reliable and efficient a-posteriori estimates in the /ip-setting. At 
the core of our analysis is the reduction of the problem on the analysis of the associated boundary value 
. problem. We start from the analysis of Wohlmuth and Melenk and combine this with our a-posteriori 

' estimation framework to obtain eigenvalue/eigenvector approximation bounds. 

O 

Q 

Accurate computation of eigenvalues and eigenvectors of differential operators defined in planar re- 
<^ , gions has attracted considerable attention recently. A central paper in this body of work is the 2005 
^ I contribution of Trefethen and Betcke [7] on computing eigenpairs for the Laplacian on a variety of planar 
domains, by the method of particular solutions. This approach produced highly accurate eigenvalues — 
correct to 13 or 14 digits in some cases — but the approach is limited in its application scope to differential 
operators whose highest order coefficients are constant and lower order coefficients are analytic, see the 
discussion from [9]. In particular this means that handling anisotropic diffusion operators is excluded. 
For further discussion of recent research in this area see [5, 6, 19]. 
j> ■ This limitation excludes many interesting eigenvalue model problems for composite materials, such as 
^ [ those which are of interest for methods of nondestructive sensing (cf. [1, 2]). Our interest in problems of 
^ ■ this sort is motivated by considerations of photonic crystals and related problems, cf. [3, 10]. Although 
O ' such problems are not directly addressed in this work, we do consider examples which have jump 
psj ■ discontinuities in the coefficients of the highest and lowest order derivatives and therefore capture some 
of the computational difficulties which arise in photonic crystal applications. In [11], we used an hp- 
adaptive discontinuous Galerkin method, with duality-based (goal-oriented) adaptive refinement, to 
efficiently produce eigenvalue approximations having at least 10 correct digits for several model problems, 
including those with discontinuous coefficients. 

Our experience thus far indicates that such hp-DG methods provide the most efficient means of 
■ computing eigenvalues in the discontinuous-coefficient case in terms of fiops-per-correct-digit. However, 
the structure of DG-methods is such that only limited results are available on the accuracy of computed 
eigenvector approximations. This is, in part, due to the difficulty in choosing an appropriate norm for 
the analysis. The analytical framework that we have developed elsewhere for lower-order continuous 
elements ([4,12]) uses native Hilbert space norms in an essential way, so standard DC norms appear 
very difficult to incorporate in this kind of analysis. 

Because it is straight-forward to apply our framework in the analysis of approximations of eigenvec- 
tors of low regularity, H^~^°' for some (small) a > 0, as well as invariant subspaces corresponding to 
degenerate eigenvalues (those having multiplicity greater than one), it seems useful to develop a contin- 
uous /ip-adaptive scheme based on this approach. The aim is that a more robust theory might soon be 
complemented with computational efficiency which is competitive with its DG counterpart. The present 
work is a first significant step in that direction. 
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The rest of this paper is organized as follows. In Section 2 we introduce our model problem and basic 
notation related to its continuous and discrete versions, as well as some basic theory related to such 
eigenvalue problems. The notion of approximation defects and their relevance is discussed in Section 3, 
with results from [4, 12] extended for use in the present context. These extensions make possible the 
incorporation of results from [16, Section 3], which pertain to boundary value problem error estimation, 
to obtain efficient and reliable estimates of eigenvalue approximations, which is discussed in Section 4. 
Also in this section we present a sin B type result for the accuracy of eigenvectors — to assess the accuracy 
of the angle operator we use the Hilbert-Schmidt norm. Section 5, which constitutes the bulk of the 
paper, is devoted to numerical experiments on a variety of different kinds of problems to assess the 
practical behavior of the proposed approach. 

2. Model Problem and Discretization 

Let 1] C be a bounded polygonal domain, possibly with re-entrant corners, and let dVLj:, C dVt 
have positive (ID) Lebesgue measure. We define the space "H = {u G H^iVL) : v\g^^ = 0}, where these 
boundary values are understood in the sense of trace. We are interested in the eigenvalue problem: 

(2.1) Find (X^tp) eRxU so that B{tp, v) = A(V^, v) and t/j ^ for all v eU . 
Here we have assumed 

(2.2) B{w,v)= / AVw -Vv + cwvdx, 

Jn 

and 

(2.3) {w,v) = / wvdx 



n 



is the standard inner-product. We will also assume that the diffusion matrix A is piecewise constant 
and uniformly positive definite a.e., the scalar c is also piecewise constant and non- negative. These as- 
sumptions are sufficient to guarantee that there are constants /3o, /3i > such that B{v, w) < ||i||w||i 
and B[v] = B{v,v) > (3o\\v\\l for all v,w El-L. In other words i?(-, ■) is an inner product on "H, whose 
induced "energy"-norm ||| ■ ||| is equivalent to || ■ ||i. The numbers /3o and (3i are referred to, respectively, 
as the coercivity and continuity constants for B. 

Here and elsewhere, we use the following standard notation for norms and seminorms: for A; G N and 
S <zVt we denote the standard norms and semi- norms on the Hilbert spaces H^{S) by 

(2.4) IK'lli5= E ll^^^ll^ K'll5=Ell^"^ll^' 

|ci|<A,' |o|=fc 

where \\ ■ \\s denotes the norm on S. Additionally, we define ||| ■ by 

(2.5) |||f ill = AVv -Vv + cv"^ dx , 

Js 

recognizing that this may be a semi-norm. When S = Q we omit it from the subscript. Our assumptions 
on A and c guarantee that there are local constants /3oS)7i5 > such that /Sosl^lis < ili'ill < Asll'^lli s, 
and the seminorm in the lower bound can be replaced with the full norm (after modifying f3os if necessary) 
if c(x) > cs > on S. 

2.1. Notational conventions for eigenvalues and eigenvectors. The variational eigenvalue prob- 
lem (2.1)-(2.3) is attained by the positive sequence of eigenvalues 

(2.6) < Ai < A2 < • ■ ■ < Ag < ■ ■ ■ 

and the sequence of eigenvectors {'ipi)ieN such that 

(2.7) B{'ijji,v) = Xi{ilJi,v), Vt; G "H, and = Sij . 



Here we have counted the eigenvalues according to their multiphcity and we will also use the notation 
ipi -L 'ipj when {ipi,ipj) = (when i 7^ j). Furthermore, the sequence (Aj)jgN has no finite accumulation 
point; and due to the Peron-Frobenius theorem we know that, in the case in which f2 is a path- wise 
connected domain, the inequality Ai < A2 holds and the eigenvector ipi can be chosen so that ipi is 
continuous and ipi > holds pointwise in fl. We will also use the notation 

Spec^ := {Xi : i e N} 

to denote the spectrum of the variational eigenvalue problem (2.7) and we use 

9Jl(A) := span{V^ : = 1, and 5(V^, 0) = A(V^, 0), V0 G 7/} 

to denote the spectral subspace associated to A G Spec^. For variational eigenvalue problems like 
(2.2) and (2.7) the subspaces 971(A), A G Spec^ are finite dimensional. Furthermore, let Ex be the 
orthogonal projection onto 9Jl(A) then 

E ^> = ^ 

and the spaces DJt{\) = Ran Ex and 9Jl(/x) = Rani^^^ are mutually orthogonal for A,/i G Spec^ and 
A ^ /i. 

Finally, note that 

5(^,0) = Yl ^(^'^A0), V^,0G'H 

AeSpcc(A) 

and so we obtain an alternative representation of the energy norm 

(2.8) \l^f = Bii;,^)= Yl K^^Ex^)^ ^ e 

AGSpcc(A) 

2.2. Discrete eigenvalue/eigenvector approximations. We discretize (2.1) using /ip-finite element 
spaces, which we now briefiy describe. Let T = Th he a triangulation of Q with the piecewise-constant 
mesh function h : Th ^ (0)1); h{K) = diam(ii') for K G Th- Throughout we implicitly assume 
that the mesh is aligned with all discontinuities of the data A and c, as well as any locations where 
the (homogeneous) boundary conditions change between Dirichlet and Neumann. Given a piecewise- 
constant distribution of polynomial degrees, p : 7h — > N, we define the space 

V = V^ = {v enn cm : G Pp(i^) for each K e %} , 

where Pj is the collection of polynomials of total degree no greater than j on a given set. Suppressing 
the mesh parameter h for convenience, we also define the set of edges S in T, and distinguish interior 
edges Si, and edges on the Neumann boundary (if there are any). Additionally, we let T(e) denote 
the one or two triangles having e G £^ as an edge, and we extend p to £^ by p{e) = maxKer{e) p{K) . As 
is standard, we assume that the family of spaces satisfy the following regularity properties on Th and p: 
There is a constant 7 > for which 
(CI) 7-^[/i(K)]2 < area(/s:) hi K eT, _ _ 

(C2) -f-\p{K) + 1) < p{K') + 1 < -fip{K) + 1) for adjacent K, K' e r,K nW ^ 0. 

It is really just a matter of notational convenience that a single constant 7 is used for all of these 
upper and lower bounds. The shape regularity assumption (CI) implies that the diameters of adjacent 
elements are comparable. 

In what follows we consider the discrete versions of (2.1): 

(2.9) Find (A, ^) G M x V" such that 5(^, v) = X{ip, v) for all v e V . 

We also assume, without further comment, that the solutions are ordered and indexed as in (2.6), with 
{ipi, ipj) = We are interested in assessing approximation errors in collections of computed eigenvalues 
and associated invariant subspaces. Let Sm = {fik}T=i ^ l*^? ^) be the set of all eigenvalues of B, counting 
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multiplicities, in the interval (a, 6), and let Sm = span{0fc}^^ be the associated invariant subspace, 
with {(f)i,(f)j) = 6ij. The discrete problem (2.9) is used to compute corresponding approximations Sm = 
{AfclfcLi and Srn = span{0fc}'^!Li, with (0^,0^) = 5ij. 

Remark 2.1. When Sm consists of the smallest m eigenvalues, we use the absolute labelling Sm = {A^}^^ 
and Sm = spanj-j/^fclfcLi instead of the relative labelling involving (/i^, (pk)', and the analogous statement 
holds for the discrete approximations Sm and Sm- This distinction is used in some of our results, such 
as Theorems 3.1 and 3.3. 



3. Approximation Defects 

3.1. Approximation defects. Let the finite element space V cTihe given and let Sm and 5*^ be the 
approximations which are computed from V. We define the corresponding approximation defects as: 

, . lll^(/)-^(/)llP 

(3.1) Vi{Sm)= ma;X mm ^771^ , 

dim5=m-i+l J^O 

where u{f) and u(f) satisfy the boundary value problems: 

(3.2) B{u{f), v) = (/, v) for every v eK 

(3.3) B{u{f), v) = (/, v) for every v E V . 

In Theorems 3.1 and 3.3 below, we state key theorems from [4, 12], which show that these approximation 
defects would yield ideal error estimates for eigenvalue and eigenvector computation if they could be 
computed. This motivates the use of a posteriori error estimation techniques for boundary value 
problems to efficiently and reliably estimate approximation defects. In [4, 12], we used hierarchical bases 
to estimate the approximation defects when V was the space of continuous, piecewise affine functions. 
In Section 4 we show how to utilize the theory of residual based estimates for /ip-finite elements from 
[16] in a similar fashion. 

The following result concerns approximations Sm and Su of the (complete) lower part of the spectrum. 
This is the reason why we have capitalized the dimension parameter M G N, which is associated to the 
cluster of lowermost eigenvalues. As opposed to a given cluster of eigenvalues contained in the interval 
(a, 6). 

Theorem 3.1. Assume that Am < Am+i, and let Sm be the span of first MEN eigenvectors of (2.9). // 



'S'm = span{^/;i, ■ ■ ■ , ?/;m} is such that ''"^ < 



then 



f M M J , M 

(3.4) ^ ^'(^m) < < Cm Y: vKS.^. 

2Am Aj .^-^ 

The constant Cm depends solely on the relative distance to the unwanted component of the spectrum (e.g. 

Am— Am+1 ) 
Am+Am+1 ' ■ 

The constant Cm is given by an explicit formula which is a reasonable practical overestimate, see 
[4, 12] for details. A similar results holds for the eigenvectors. We point the interested reader to [12, 
Theorem 4.1 and equation (3.10)] and [4, Theorem 3.10]. 

Remark 3.2. Although Ai < A2 for the particular problems we consider numerically in the present 
work, much of the theory carries over to problems where Q is not pathwise connected, or the boundary 
conditions are periodic (as examples). In these cases the Peron-Frobenius theorem does not apply, and 
it is quite possible that the smallest eigenvalue is degenerate. If this is the case, and Ai = Am, then the 
constant Ai/2Am in (3.4) can be replaced by 1. 
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An important feature of these ideal estimates is that they are asymptotically exact, both as eigenvector 
as well as as eigenvalue estimators, as the following theorem indicates in the case of a single degenerate 
eigenvalue and its corresponding invariant subspace. 



Theorem 3.3. Let Xq be a degenerate eigenvalue of multiplicty m, Ag_i < = Xq+m-i < ^q+m- 
Let Sm = Sm{T) = span(0fc) G V = V{T) be the computed approximation of the invariant subspace 
corresponding to Xg. Then, taking the pairing of eigenvectors 0, and Ritz vectors (pi as in [12], we have 



Em Ifii—Xql \-^m B[t 

1=1 hi ^ Z^i=l 



(3.5) lim ^1 = 1 , lim ^'=' = 1 . 

3.2. A relationship with the residual estimates for a Ritz vector basis. This section addresses 
the issue of the computability of r]i{Sm) by relating these quantities to the standard dual energy norm 
estimates of the residuals associated to the Ritz vector basis of Sm- 

In our notation the energy norm was denoted by ||| ■ ||| and we use u{-) and u(-) to denote the solution 
operators from (3.2) and (3.3). We assume are the Ritz vectors from Sm., then for i,j = 

1, . . . , m, we define the matrices 

(3.6) Eij = B{u{(j)i) - u{4>i), u{(j)j) - u{4>j)) 

(3.7) Gij = B{u{^i),u{^,)). 

These matrices were introduced in [12] under the name of the error and the gradient matrix. It was 
shown in [12] that rji^Sm) = Xi{E,G), where Xi{E,G) < ■■■ < Xm{E,G) are the eigenvalues of the 
generalized eigenproblem for the matrix pair {E, G). 

We further assume that 0^, i = 1, ■ ■ ■ , m are among the Ritz vectors from the finite element subspace 
V, V D Sm from (3.3). The identity (3.6) implies that E is a. Gram matrix for the set of vectors 
u{(j)i) — u{4>i), i = 1, . . . ,m. If we assume that Sm does not contain any eigenvectors, then we conclude 
that E must be positive definite matrix. Furthermore, it always holds 

r^^iSm) = X,{G-'/'EG-'/') 

(3.8) Eii = ^^'^\\u{fii(i)i) - u{fii(i)i)f , i = l,...,m 

D,,<G<{1 + ^i)D^, 

where D^^ = diag{fi]^^, . . . , fi^) and = ||-D^ ""^^^(G — Z^^)-D^ ^^^||. Let us note that Di is the relative 
estimate, so it is expected that even for very crude finite element spaces V we have < 1. 
Now compute 

m m 

Y.X,{D-'I^ED~''^) = tr{D~'l-ED~''') = J^E.f^^, 
1=1 1=1 

and so conclude that 

^ m m. m 

(3.9) — — E^^fi^ < vKS,n) < E,,fli . 

1 + '-r 

1=1 1=1 1=1 

We summarize this considerations — using the identity (3.8) — in the following lemma. 
Lemma 3.4. It holds that 

^ m mm 

^^'^^"^ 1 + S) ^f^i^t^ii^i^i) ~ Hi^i<Pi)f < ^^ViiSm) < ^ /i^^ |||u(/ii0i) - u{fli4>i)f ■ 

1=1 1=1 i=l 
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4. /ip-ERROR Estimation and Adaptivity in the Eigenvalue Context 



Using Lemma 3.4, we have reduced the problem of estimating the approximation defects, and hence 
the error in our eigenvalue/eigenvector computations, to that of estimating error in associated boundary 
value problems. In particular, we must estimate — for each Ritz vector, where = 

span{0i, . . . , (pm} is our approximation of Sm = span{0i, . . . , (pm}- We modify key results from [16], 
which were stated only for the Laplacian, to our context. The identity u{fii(j)i) = (pi, makes our job 
easier. We define the element residuals Ri for K G T, and the edge (jump) residuals rj for e E S,hy 



(4.1) Ri\K = fii^i- c^i + V ■ AV^i , 

-{AV^i)l^■nK-{AV(p^)\^,■nK' , e G 
-(AV0i)l^ UK , e e£N 



(4.2) 



For interior edges e & Sj, K and K' are the two adjacent elements, having outward unit normals uk and 
iiK', respectively; and for Neumann boundary edges e E £n (if there are any), K is the single adjacent 
element, having outward unit normal n^- We note that i? is a polynomial of degree no greater than 
p{K) on K, and r is a polynomial of degree no greater than p{e) on e. 

Our estimate of ef = YliKeT^'i^^) ~ ~ "^(Ai^OIlP is computed from local quantities, 

where £i{K) and £n{K) denote the interior edges and Neumann boundary edges of K, respectively. 
An inspection the proof of [16, Lemma 3.1] (which was stated for the Laplacian) makes the following 
assertion clear. 

Lemma 4.1. There is a constant C > depending only on the hp-constant 7 and the coercivity constant 
(3o, such that |||u(/ij0j) — u(/ij</)j)|||^ < Cej. 

A few remarks are in order concerning the lemma above and how it relates to [16, Lemma 3.1]. First, 
the bound in [16, Lemma 3.1] includes an additional term involving the difference between the righthand 
side (in our case (pi) and its projection on K into a space of polynomials. This additional term only 
arises in their result because they have chosen to use the projection of the righthand side, instead of the 
righthand side itself, to define the element residual (here called Ri). They do this in order to employ 
certain polynomial inverse estimates, which hold in our case outright because our righthand sides are 
piecewise polynomial. Their result also involves a parameter a G [0, 1], which we have taken to be 0. 
The result [16, Lemma 3.1] is based on Scott-Zhang type quasi-interpolation, which naturally gives rise 
to errors measured in H^. Mimicking their arguments with our indicator, one would arrive at a result 
of the form 



u{fii(pi) - u{iii(f)i)l < Cei\\u{fLi(pi) - u{iii(pi)\\ 



where C depends only on 7. The constant in the coercivity bound /Soll'^lli ^ il'^li^ enters Lemma 4.1 at 
this final stage. Similarly, a careful reading of the proofs of [16, Lemma 3.4 and 3.5] show that their 
efficiency results are readily extended to elliptic operators of the type considered here. 

Lemma 4.2. For any e > 0, there is a constant c = c(e) > depending only on the hp-constant 7 and 
the global continuity constant f3i, such that eji^K) < cp']^'^'^\u[{ii(pi) — u{jli(pi)\'^^. 

Here, ujk is the patch of elements which share an edge with K. The global continuity constant /3i could 
be replaced in Lemma 4.2 by a local continuity constant /Siojj^ if desired. 

Remark 4.3. The p-dependence in local efficiency bound of Lemma 4.2 is unfortunately unavoidable in 
the proof, and would suggest decreased efficiency of the estimator as pk is increased if this estimate 
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were sharp. Our numerical experiments do seem to indicate that efficiency does, in fact, decrease under 
/ip- refinement, but that this decrease is modest in practical computations. 

With these results we now state the main theorem. 

Theorem 4.4. Under the assumptions of Theorem 3.1, we have the following upper- and lower-hounds 
on eigenvalue error, 

M M J , M 

(4.4) Y: k'e^ < ^ E ■ 

i=l i=l i=l 

The constant Ci depends solely on the ratio Ai/(2A2), the hp -regularity constant 7, the continuity con- 
stant Pi, and the maximal polynomial degree p = maxK^rPiK). The constant C2 depends solely on 
the relative distance to the unwanted component of the spectrum, the hp-regularity constant 7 and the 
coercivity constant Pq. 

Proof. These assertions follow directly from Theorem 3.1, Lemma 3.4, and Lemmas 4.1 and 4.2. Q.E.D. 

Remark 4.5. It is relative local indicators fi~^e^{K) which will be used to mark elements for refinement, 
as will be described in Section 5. 

A similar result holds for the eigenvectors and eigenspaces. Let now 

A<Am, AGSpec^ 

be the orthogonal projection onto the eigenspace belonging to the first M eigenvalues of the form B as 
given in Theorem 4.4. Let also || ■ \\s2 be the Hilbert-Schmidt norm on the ideal of all Hilbert-Schmidt 
operators, see [18]. We now have the eigenvector result. 

Theorem 4.6. Under the assumptions of Theorem 3.1, we have the following upper- and lower-hounds 
on eigenfunction error. 



(4.5) ||sine(E(AM),^M)||& <CM,r 



i=l 



The constant Cm,t depends solely on the relative distance to the unwanted component of the spectrum 
(e.g. ^^^^Jj the hp-regularity constant'-]' and the continuity constant f3i. 

5. Experiments 



In the numerical experiments we illustrate the efficiency of the estimator (4.4) on several problems of 
the general form 

(5.1) Li/j = Xi/j inn , ii^:-!! = i , 

for a second-order, linear elliptic operator C, where homogeneous Dirichlet or Neumann conditions are 
imposed on the boundary. Plots are given of the total relative error, its a posteriori estimate, and the 
associated effectivity quotient, shown, respectively, below: 



These are plotted against the square- root of the size of the discrete problem DOFs = dim(V"^). For 
most of our problems, the exact eigenvalues are unknown, so we take highly accurate computations on 
very large problems as "exact" for these comparisons. 

In all simulations we used an /ip-adaptive algorithm in order to get the best convergence possible. To 
drive the /ip-adaptivity we use the element-wise contributions to the quantity Y^i=i •^r^^o to provide 



Figure 1. Some of the domains used for the experiments. 



local error indicators. Then, we apply a simple fixed- fraction strategy to mark the elements to adapt. 
For each marked element, the choice of whether to locally refine it or vary its approximation order is 
made by estimating the local analyticity of the computed eigenf unctions in the interior of the element 
by computing the coefficients of the L^-orthogonal polynomial expansion (cf. [8]). 

5.1. Dirichlet Laplacian on the Unit Triangle. As a simple problem for which the eigenvalues and 
eigenfunctions are explicitly known (cf. [15]), we consider the problem where: L = —A, Q is equilateral 
triangle of having unit edge-length, and = on dQ. The eigenvalues can be indexed as 

>^mn = —f—V^ +mn + n ) , 

y 

and we refer interested readers to [15] for explicit descriptions of the eigenfunctions. 




Figure 2. Errors and error estimates. Triangle problem. 

In Figure 2 we plot the total relative error for the first four eigenvalues, together with the associated 
error estimate; and in Figure 3 we plot the effectivity quotient. It is clear that the convergence is 
exponential in this case, and that the effectivity undergoes a mild degradation as the problem size 
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0.042 




5 10 15 20 25 30 

Figure 3. Effectivity index. Triangle problem. 

increases. This modest decrease in effectivity is in line with Remark 4.3, and it is also seen in our 
remaining experiments. 




Figure 4. Errors and error estimates. Triangle with hole. 

5.2. Dirichlet Laplacian on the Unit Triangle with on a Hole. We now consider the problem 
where L = —A, Q is the equilateral triangle having edge-length 2 with an equilateral triangle having 
edge-length 1/2 removed from its center (see Figure 1), and ip = on dQ. For such a problem, it is 
expected that some of the eigenfunctions will have an r^/^-type singularity at each of the three interior 
corners, where r is the distance to the nearest corner. In this case, the exact eigenvalues are unknown, 
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Figure 5. Effectivity index. Triangle with hole. 



so we computed the following reference values of them on a very large problem: 40.4650426 for the first 
eigenvalue and 43.4868466 for the second and third, which form a double eigenvalue. These values are 
accurate at least up to le-6. 

In Figure 4 we plot the relative error and error estimates together, for the first three eigenvalues, 
and in Figures 5 we plot the corresponding values of the effectivity quotient. We again see exponential 
convergence and a modest deterioration of effectivity. 

5.3. Square Domain with Discontinuous Reaction Term. For this pair of problems we take Q = 
(0, 1)^, Vip ■ n = on dQ, and Cip = —/S.ip + hVmd ■ i^, where Vmd is the characteristic function 
of the touching squares labelled A^i in Figure 6. We consider two values of the constant parameter, 
K = 10, 100. It is straightforward to see that the corresponding bilinear form is an inner-product in this 
case (no zero eigenvalues), and that all eigenf unctions are at least in H^. 



M2 






M2 



Figure 6. A modification of the touching squares example of M. Dauge. 
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For K = 10, we have in Figure 7 the total relative error and error estimates for the first four eigenvalues; 
and the effectivity quotient is given in Figure 8. For these simulations we used the following reference 
values for the first four eigenvalues, which are le-8 accurate: 4.150242455, 10.706070962, 18.779725462, 
25.150325247. The analogous plots for the first four eigenvalues in the case k, = 100 are given in Figure 9 
and Figure 10. For these simulations, we used the following reference values for the first four eigenvalues, 
which are le-8 accurate: 13.210576406, 13.990033964, 60.294151672, 64.840268299. In both cases we 
see apparent exponential convergence, and reasonable effectivity behavior. It seem clear from the error 
plots that for both values of k the convergence is exponential. 




Figure 7. Errors and error estimates. Discontinuous reaction term, k = 10. 

5.4. Square Domain with Discontinuous Diffusion Term. Using the domain Q = (0, 1)^, parti- 
tioned into regions A^i and Ai^ as in Figure 6, and homogeneous Dirichlet conditions = on dfl, we 
consider the operator C = —V ■ (aV), where a = 1 in AI2 and a in A^i may vary. Such problems can 
have arbitrarily bad singularities at the cross-point of the domain depending on the relative sizes of a 
in the two subdomains — see, for example, [13, 14] and [17, Example 5.3]. 

We have considered two values for a in Aii. 10 and 100. Since the exact eigenvalues are not avail- 
able, we computed the following three reference values for the first three eigenvalues when a = 10: 
64.226529416, 75.028156269, 141.161506328; and the following three reference values for the first three 
eigenvalues when a = 100: 77.800981966, 78.564198245, 193.916538067. All reference values are at least 
le-8 accurate. The relative error and effectivity plots for both cases are given in Figures 11-14, and 
again we see apparent exponential convergence. 

5.5. Square Domain with a Slit. For this problem, C = —A and Q = (0, 1)^ \ 5", where 5* = 
{(x, 1/2) : 1/2 < X < 1}; this is pictured in Figure 1, with S as the dashed segment. Homogeneous 
Neumann conditions are imposed on both "sides" of S and homogeneous Dirichlet boundary conditions 
are imposed on the rest of the boundary of Q. For this example we used the following reference values for 
the first four eigenvalues with in brackets the corresponding accuracy: 20.739208802(8), 34.485320(5), 
50.348022005(8), 67.581165196(8). 

To give some indication of the nature of the eigenfunctions in the interior, we briefly consider a related 
problem where fl is the unit disk with a slit along the positive x-axis, as pictured in Figure 1, with 
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Figure 8. Effectivity index. Discontinuous reaction term, k, = 10. 
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Figure 9. Errors and error estimates. Discontinuous reaction term, k = 100. 



the same boundary conditions. In this case, the eigenvalues and eigenf unctions are known exphcitly. 
For k > and m > 1, let be the m*^ positive root of the first-kind Bessel function Jk/2- It is 
straightforward to verify that, up to renormalization of eigenf unctions, the eigenpairs can be indexed by 



A 



km 



ipkm = Jk/2{zkmr)sm{k9/2) 



k,m eN 



1 /2 

We see that tpkm ~ sin(A;6'/2) (^^f^) as r — 0, so singularities of types r'^/^ occur infinitely many times 
in the spectrum. The strongest of these singularities is of type r^/^, and it occurs in the eigenfunction 
associated with the second eigenvalue, for example. The same asymptotic behavior of the eigenfunctions 
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Figure 10. Effectivity index. Discontinuous reaction term, k, = 100. 




20 40 60 80 100 



Figure 11. Errors and error estimates. Discontinuous diffusion term, k = 10. 

near the crack tip is expected for the square and circular domains, and in Figure 15 we show contour plots 
of the second eigenvalue in both cases. For the circular domain, the second eigenvalue and function can, in 
principle be obtained to arbitrary precision using a computer algebra system. Using Mathematica, we 
computed the second eigenvalue (the smallest positive root of J1/2) to 20 digits, 9.869604401089358619, 
and generated the corresponding contour plot of the eigenf unction. 

In Figure 16 we plot the total relative errors and error estimates for the first four eigenvalues, and in 
Figure 17 the individual eigenvalue errors are shown. It is clear from the second of these figures that the 
second, which corresponds to the most singular eigenfunction, clearly has the worst convergence rate 
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Figure 12. Effectivity index. Discontinuous diffusion term, k = 10. 
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Figure 13. Errors and error estimates. Discontinuous diffusion term, k = 100. 



(as expected), and that this is what "spoils" the convergence of the cluster of the first four eigenvalues. 
This becomes even more apparent when Figure 18, which corresponds to the second eigenvalue alone, 
is compared with Figure 16 — they are nearly identical. 
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Figure 14. Effectivity index. Discontinuous diffusion term, n = 100. 





Figure 15. Second eigenfunction for the slit circle (top) and slit square. 
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